options(warn=2)
options(error=recover)
Sys.setenv(TZ='EST')
options(show.error.locations=TRUE)

dir <- getwd()

cat('original.path >> ')
orig.data.path <- readLines(n=1)
setwd(orig.data.path)
orig.data <- sapply(list.files(pattern='\\.csv$'), read.csv, header=TRUE, USE.NAMES=TRUE, simplify=FALSE)
seg.nums1 <- sapply(regmatches(names(orig.data), regexpr('\\.csv$', names(orig.data)), invert=TRUE), paste, collapse='')

cat('corrected.path >> ')
corrected.data.path <- readLines(n=1)
setwd(corrected.data.path)
corrected.data <- sapply(list.files(pattern='\\.csv$'), read.csv, header=TRUE, USE.NAMES=TRUE, simplify=FALSE)
seg.nums2 <- sapply(regmatches(names(corrected.data), regexpr('\\.csv$', names(corrected.data)), invert=TRUE), paste, collapse='')

cat('basevol path >> ')
basevol.data.path <- readLines(n=1)
setwd(basevol.data.path)
basevol.data <- sapply(list.files(pattern='\\.csv$'), read.csv, header=TRUE, USE.NAMES=TRUE, simplify=FALSE)
seg.nums3 <- sapply(regmatches(names(basevol.data), regexpr('\\.csv$', names(basevol.data)), invert=TRUE), paste, collapse='')

setwd(dir)

cat('output dir >> ')
out.dir <- file.path(readLines(n=1), 'compare_plots_vol')

if (! isTRUE(file.info(out.dir)$isdir)) {
    dir.create(out.dir)
} else {
    unlink(out.dir, recursive=TRUE) 
}

seg.nums <- if (length(seg.nums1) > length(seg.nums2)) seg.nums2 else seg.nums1

cat('enter time boundaries? >> ')
yes.no <- readLines(n=1)
if (grepl('y', yes.no)) {
    cat('enter report start >> ')
    begin.time <- as.POSIXct(readLines(n=1))
    cat('enter timestep >> ')
    timestep <- as.difftime(as.numeric(readLines(n=1)), units='secs')
    cat('enter event start time >> ') 
    start.dtime <- as.POSIXct(readLines(n=1))
    cat('enter event end time >> ')
    end.dtime <- as.POSIXct(readLines(n=1))
   
    timesteps <- seq(from=begin.time, to=end.dtime, by=timestep)
    subset.indices <- which(!(timesteps < start.dtime))
} else {
    subset.indices <- sequence(nrow(orig.data[[1]]))
}


plotter <- function(orig.data, corrected.data, basevol.data, seg.num, out.dir) {
    ylim <- c(min(orig.data$vol, corrected.data$vol), max(orig.data$vol, corrected.data$vol))
    with(orig.data, {
        main <- paste(seg.num, 'original vs. corrected volume', sep=' - ')
        basevol <- mean(vol[1:721], na.rm=TRUE)
        png(filename=file.path(out.dir, paste(main, 'png', sep='.')), height=1000, width=1000, pointsize=12)
        par(mar=c(5,6,4,6) + 0.1)
        yvals <- vol[subset.indices]
        #ylim <- ylim=c(min(yvals, basevol), max(yvals, basevol))
        plot(timestep[subset.indices], yvals, main=main, type='b', 
             lwd=2, ylab='volume (m3)', yaxt='n', cex=.2, ylim=ylim)
        axis(side=2, at=pretty(vol, n=25))
        abline(h=basevol, col='blue')
    })

    with(corrected.data, {
        lines(timestep[subset.indices], vol[subset.indices], type='b', lwd=1, col='red', cex=.2)
    })

    #with(basevol.data, {
    #    lines(timestep[subset.indices], vol[subset.indices], type='b', lwd=1, col='green', cex=.2)
    #})

    legend('topleft', legend=c('original volume', 'corrected volume'), col=c('black', 'red'), lty=1)
    graphics.off()
}

#pdf(file=file.path(out.dir, 'out.pdf'))
mapply(plotter, orig.data[1:length(seg.nums)], corrected.data[1:length(seg.nums)], basevol.data[1:length(seg.nums)], 
       seg.nums, out.dir)
#graphics.off()

